!**********************************************************
!** This file contains the UBC-IO subroutines            **
!**********************************************************
!* Created by Zhang Sheng at CUGB                        **
!*      version: 2015.5.20                               **
!*      email:zhangsheng@cugb.edu.cn                     **
!**********************************************************
!  Need include 'datio.f90'                              **
!**********************************************************
!>>>>>>>>>           The input part            <<<<<<<<<<<<
!**********************************************************
!   --------------------  common ------------------------
! (  1 ) read_mesh_size(meshfile,ne,nn,nv)
! (  2 ) read_mesh_dat(meshfile,ne,nn,nv,e0,n0,v0,de,dn,dv)
! (  3 ) read_topo_npt(topofile,npt)
! (  4 ) read_topo_dat(topofile,npt,ee,nn,ele)
! (  9 ) read_phys_dat(denfile,ne,nn,nv,den)
! ( 10 ) read_bound_dat(denfile,ne,nn,nv,lden,uden)
! ( 11 ) read_weight_dat(denfile,ne,nn,nv,ws,we,wn,wz)
!
!   --------------------  gravity  ------------------------
! (  5 ) read_gloc_npt(gobsfile,npt)
! (  6 ) read_gloc_dat(gobsfile,npt,oee,onn,oele)
! (  7 ) read_gobs_npt(gobsfile,npt)
! (  8 ) read_gobs_dat(gobsfile,npt,gee,gnn,gele,grv,gerr)
!
!   --------------------  magnetic ------------------------
! ( 12 ) read_mloc_mag(mlocfile,incl,decl,geomag,aincl,adecl,idir)
! ( 13 ) read_mloc_npt(mlocfile,npt)
! ( 14 ) read_mloc_dat(mlocfile,npt,oee,onn,oele,ai,ad)
! ( 15 ) read_mobs_mag(mobsfile,incl,decl,geomag,aincl,adecl,idir)
! ( 16 ) read_mobs_npt(mobsfile,npt)
! ( 17 ) read_mobs_dat(mobsfile,npt,oee,onn,oele,ai,ad)
! 
!**********************************************************
!>>>>>>>>>           The output part           <<<<<<<<<<<<
!**********************************************************
!   --------------------  common ------------------------
! ( 18 ) write_mesh_sub(meshfile,ne,nn,nv,e0,n0,v0,de,dn,dv) 
! ( 19 ) write_topo_sub(topofile,npt,ee,nn,ele) 
! ( 20 ) write_phys_sub(denfile,ne,nn,nv,den) 
! ( 21 ) write_bound_sub(denfile,ne,nn,nv,lden,uden) 
! ( 22 ) write_weight_sub(denfile,ne,nn,nv,ws,we,wn,wz) 
!
!   --------------------  gravity  ------------------------
! ( 23 ) write_gloc_sub(glocfile,npt,oee,onn,oele)
! ( 24 ) write_gobs_sub(gobsfile,npt,gee,gnn,gele,grv,gerr)
!
!   --------------------  magnetic ------------------------
! ( 26 ) write_mloc_sub(mlocfile,incl,decl,geomag,aincl,adecl,idir,npt,oee,onn,oele,ai,ad)
! ( 27 ) write_mobs_sub(mobsfile,incl,decl,geomag,aincl,adecl,idir,npt,oee,onn,oele,ai,ad,mag,merr)
! 
!**********************************************************
!///////////////////// input part /////////////////////////
!**********************************************************
!(1) input m,n from infile
subroutine read_mesh_size(meshfile,ne,nn,nv)
character*(*) meshfile
character*200 tmp
integer ne,nn,nv
integer nut
call open_or_not(nut)
open(nut,file=meshfile,status='old')
do
  read(nut,'(a200)') tmp
  tmp=adjustl(tmp)
  if(tmp(1:1).ne.'!') then
    read(tmp,*) ne,nn,nv
    exit
  endif
enddo
close(nut)
end subroutine
!**********************************************************
!(2) input e0,n0,v0,de,dn,dv from meshfile
subroutine read_mesh_dat(meshfile,ne,nn,nv,e0,n0,v0,de,dn,dv)
character*(*) meshfile
character*200 tmp
integer ne,nn,nv
real e0,n0,v0
real de(ne),dn(nn),dv(nv)
integer nut
call open_or_not(nut)
open(nut,file=meshfile,status='old')
do
  read(nut,'(a200)') tmp
  tmp=adjustl(tmp)
  if(tmp(1:1).ne.'!') then
    exit
  endif
enddo
read(nut,*) e0,n0,v0
read(nut,*) (de(i),i=1,ne)
read(nut,*) (dn(i),i=1,nn)
read(nut,*) (dv(i),i=1,nv)
close(nut)
end subroutine
!**********************************************************
!(3)input number of point from topo file
subroutine read_topo_npt(topofile,npt)
character*(*) topofile
character*200 tmp
integer npt
call open_or_not(nut)
open(nut,file=topofile,status='old')
do
  read(nut,'(a200)') tmp
  tmp=adjustl(tmp)
  if(tmp(1:1).ne.'!') then
    read(tmp,*) npt
    exit
  endif
enddo

close(nut)
end subroutine
!**********************************************************
!(4)input ele dat from topo file
subroutine read_topo_dat(topofile,npt,ee,nn,ele)
character*(*) topofile
character*200 tmp
integer npt
integer i
real ee(npt),nn(npt),ele(npt)
integer nut
call open_or_not(nut)
open(nut,file=topofile,status='old')
do
  read(nut,'(a200)') tmp
  tmp=adjustl(tmp)
  if(tmp(1:1).ne.'!') then
    exit
  endif
enddo
do i=1,npt
  read(nut,*) ee(i),nn(i),ele(i)
enddo
close(nut)
end subroutine
!**********************************************************
!(5)input number of point from obs_loc file
subroutine read_gloc_npt(obsfile,npt)
character*(*) obsfile
integer npt
call read_topo_npt(obsfile,npt)
endsubroutine
!**********************************************************
!(6)input obs dat from obs_loc file
subroutine read_gloc_dat(obsfile,npt,oee,onn,oele)
character*(*) obsfile
integer npt
real oee(npt),onn(npt),oele(npt)
call read_topo_dat(obsfile,npt,oee,onn,oele)
end subroutine
!**********************************************************
!(7)input grv number of point from obs_grv file
subroutine read_gobs_npt(gobsfile,npt)
character*(*) gobsfile
integer npt
call read_topo_npt(gobsfile,npt)
endsubroutine
!**********************************************************
!(8)input grv obs dat from obs.grv file
subroutine read_gobs_dat(gobsfile,npt,gee,gnn,gele,grv,gerr)
character*(*) gobsfile
character*200 tmp
integer npt
real gee(npt),gnn(npt),gele(npt),grv(npt),gerr(npt)
integer nut
call open_or_not(nut)
open(nut,file=gobsfile,status='old')
do
  read(nut,'(a200)') tmp
  tmp=adjustl(tmp)
  if(tmp(1:1).ne.'!') then
    exit
  endif
enddo
do i=1,npt
  read(nut,*) gee(i),gnn(i),gele(i),grv(i),gerr(i)
enddo
close(nut)
end subroutine
!**********************************************************
!(9)input physical * dat from density or susceptibility file 
subroutine read_phys_dat(denfile,ne,nn,nv,den)
character*(*) denfile
integer ne,nn,nv
real den(ne*nn*nv)
integer nut
call open_or_not(nut)
open(nut,file=denfile,status='old')
do i=1,ne
do j=1,nn
do k=1,nv
  ijk=(i-1)*nn*nv+(j-1)*nv+k
  read(nut,*) den(ijk)
enddo
enddo
enddo
close(nut)
end subroutine
!**********************************************************
!(10)input lowden and upden dat from bounds.dat file 
subroutine read_bound_dat(boundfile,ne,nn,nv,lden,uden)
character*(*) boundfile
integer ne,nn,nv
real lden(ne*nn*nv),uden(ne*nn*nv)
integer nut
call open_or_not(nut)
open(nut,file=boundfile,status='old')
do i=1,ne
do j=1,nn
do k=1,nv
  ijk=(i-1)*nn*nv+(j-1)*nv+k
  read(nut,*) lden(ijk),uden(ijk)
enddo
enddo
enddo
close(nut)
end subroutine
!**********************************************************
!(11)input weight fro w.dat file 
subroutine read_weight_dat(wfile,ne,nn,nv,ws,we,wn,wz)
character*(*) wfile
integer ne,nn,nv
real ws(ne*nn*nv),we(ne*nn*nv),wn(ne*nn*nv),wz(ne*nn*nv)
integer nut
call open_or_not(nut)
open(nut,file=wfile,status='old')

  read(nut,*) (((ws((i-1)*nn*nv+(j-1)*nv+k),k=1,nv),j=1,nn),i=1,ne)
  read(nut,*) (((we((i-1)*nn*nv+(j-1)*nv+k),k=1,nv),j=1,nn),i=1,ne)
  read(nut,*) (((wn((i-1)*nn*nv+(j-1)*nv+k),k=1,nv),j=1,nn),i=1,ne)
  read(nut,*) (((wz((i-1)*nn*nv+(j-1)*nv+k),k=1,nv),j=1,nn),i=1,ne)

close(nut)
end subroutine
!**********************************************************
! ( 12 ) 
subroutine read_mloc_mag(mobsfile,incl,decl,geomag,aincl,adecl,idir)
character*(*) mobsfile
character*200 tmp
real incl,decl,geomag,aincl,adecl
integer idir 
integer nut
call open_or_not(nut)
open(nut,file=mobsfile,status='old')
do
  read(nut,'(a200)') tmp
  tmp=adjustl(tmp)
  if(tmp(1:1).ne.'!') then
    read(tmp,*) incl,decl,geomag
    read(nut,*) aincl,adecl,idir
    exit
  endif
enddo

close(nut)
endsubroutine

!**********************************************************
! ( 13 ) 
subroutine read_mloc_npt(mobsfile,npt)
character*(*) mobsfile
character*200 tmp
integer nut
real incl,decl,geomag,aincl,adecl
integer idir 
integer npt
call open_or_not(nut)
open(nut,file=mobsfile,status='old')
do
  read(nut,'(a200)') tmp
  tmp=adjustl(tmp)
  if(tmp(1:1).ne.'!') then
    read(tmp,*) incl,decl,geomag
    read(nut,*) aincl,adecl,idir
    read(nut,*) npt
    exit
  endif
enddo

close(nut)
end subroutine
!**********************************************************
! ( 14 ) 
subroutine read_mloc_dat(mobsfile,npt,oee,onn,oele,ai,ad)
character*(*) mobsfile
character*200 tmp
integer nut
real incl,decl,geomag,aincl,adecl
integer idir 
integer npt
real oee(npt),onn(npt),oele(npt),ai(npt),ad(npt)
call open_or_not(nut)
open(nut,file=mobsfile,status='old')
do
  read(nut,'(a200)') tmp
  tmp=adjustl(tmp)
  if(tmp(1:1).ne.'!') then
    read(tmp,*) incl,decl,geomag
    read(nut,*) aincl,adecl,idir
    read(nut,*) npt
    exit
  endif
enddo
if(idir.eq.1) then
  do i=1,npt
    read(nut,*) oee(i),onn(i),oele(i)
  enddo
elseif(idir.eq.0)then
  do i=1,npt
    read(nut,*) oee(i),onn(i),oele(i),ai(i),ad(i)
  enddo
endif  
close(nut)
end subroutine
!**********************************************************
! ( 15 ) 
subroutine read_mobs_mag(mobsfile,incl,decl,geomag,aincl,adecl,idir)
character*(*) mobsfile
real incl,decl,geomag,aincl,adecl
integer idir 
call read_mloc_mag(mobsfile,incl,decl,geomag,aincl,adecl,idir)
end subroutine
!**********************************************************
! ( 16 ) 
subroutine read_mobs_npt(mobsfile,npt)
character*(*) mobsfile
integer npt
call read_mloc_npt(mobsfile,npt)
end subroutine
!**********************************************************
! ( 17 ) 
subroutine read_mobs_dat(mobsfile,npt,oee,onn,oele,ai,ad,mag,merr)
character*(*) mobsfile
character*200 tmp
integer npt
real oee(npt),onn(npt),oele(npt),ai(npt),ad(npt),mag(npt),merr(npt)
call open_or_not(nut)
open(nut,file=mobsfile,status='old')
do
  read(nut,'(a200)') tmp
  tmp=adjustl(tmp)
  if(tmp(1:1).ne.'!') then
    read(tmp,*) incl,decl,geomag
    read(nut,*) aincl,adecl,idir
    read(nut,*) npt
    exit
  endif
enddo
if(idir.eq.1) then
  do i=1,npt
    read(nut,*) oee(i),onn(i),oele(i),mag(i),merr(i)
  enddo
elseif(idir.eq.0)then
  do i=1,npt
    read(nut,*) oee(i),onn(i),oele(i),ai(i),ad(i),mag(i),merr(i)
  enddo
endif  
close(nut)
end subroutine
!**********************************************************
!///////////////////// output part ////////////////////////
!**********************************************************
! ( 18 ) 
subroutine write_mesh_sub(meshfile,ne,nn,nv,e0,n0,v0,de,dn,dv) 
character*(*) meshfile
integer ne,nn,nv
real e0,n0,v0
real de(ne),dn(nn),dv(nv)
integer nut
call open_or_not(nut)
open(nut,file=meshfile)
write(nut,*) ne,nn,nv
write(nut,*) e0,n0,v0
write(nut,*) (de(i),i=1,ne)
write(nut,*) (dn(i),i=1,nn)
write(nut,*) (dv(i),i=1,nv)
close(nut)
end subroutine
!**********************************************************
! ( 19 ) 
subroutine write_topo_sub(topofile,npt,ee,nn,ele) 
character*(*) topofile
integer npt
real ee(npt),nn(npt),ele(npt)
integer i,nut
call open_or_not(nut)
open(nut,file=topofile)
write(nut,*) npt
do i=1,npt
write(nut,*)  ee(i),nn(i),ele(i)
enddo
close(nut)
end subroutine
!**********************************************************
! ( 20 ) 
subroutine write_phys_sub(denfile,ne,nn,nv,den) 
character*(*) denfile
integer ne,nn,nv
real den(ne*nn*nv)
integer i,j,k,ijk,nut
call open_or_not(nut)
open(nut,file=denfile)
do i=1,ne
do j=1,nn
do k=1,nv
  ijk=(i-1)*nn*nv+(j-1)*nv+k
  write(nut,*) den(ijk)
enddo
enddo
enddo
close(nut)
end subroutine
!**********************************************************
! ( 21 ) 
subroutine write_bound_sub(denfile,ne,nn,nv,lden,uden) 
character*(*) denfile
integer ne,nn,nv
real lden(ne*nn*nv),uden(ne*nn*nv)
integer i,j,k,ijk,nut
call open_or_not(nut)
open(nut,file=denfile)
do i=1,ne
do j=1,nn
do k=1,nv
  ijk=(i-1)*nn*nv+(j-1)*nv+k
  write(nut,*) lden(ijk),uden(ijk)
enddo
enddo
enddo
close(nut)
end subroutine
!**********************************************************
! ( 22 ) 
subroutine write_weight_sub(denfile,ne,nn,nv,ws,we,wn,wz) 
character*(*) denfile
integer ne,nn,nv
real ws(ne*nn*nv),we(ne*nn*nv),wn(ne*nn*nv),wz(ne*nn*nv)
integer i,j,k,ijk,nut
call open_or_not(nut)
open(nut,file=denfile)
  write(nut,'(f\)',advance='no') (((ws((i-1)*nn*nv+(j-1)*nv+k),k=1,nv),j=1,nn),i=1,ne)
  write(nut,'(f\)',advance='no') (((we((i-1)*nn*nv+(j-1)*nv+k),k=1,nv),j=1,nn),i=1,ne)
  write(nut,'(f\)',advance='no') (((wn((i-1)*nn*nv+(j-1)*nv+k),k=1,nv),j=1,nn),i=1,ne)
  write(nut,'(f\)',advance='no') (((wz((i-1)*nn*nv+(j-1)*nv+k),k=1,nv),j=1,nn),i=1,ne)
close(nut)
end subroutine
!**********************************************************
! ( 23 ) 
subroutine write_gloc_sub(glocfile,npt,oee,onn,oele)
character*(*) glocfile
integer npt
real oee(npt),onn(npt),oele(npt)
call write_topo_sub(glocfile,npt,oee,onn,oele) 
end subroutine
!**********************************************************
! ( 24 ) 
subroutine write_gobs_sub(gobsfile,npt,gee,gnn,gele,grv,gerr)
character*(*) gobsfile
integer npt
real gee(npt),gnn(npt),gele(npt),grv(npt),gerr(npt)
integer nut
call open_or_not(nut)
open(nut,file=gobsfile)
  write(nut,*) npt
do i=1,npt
  write(nut,'(5f)')  gee(i),gnn(i),gele(i),grv(i),gerr(i)
enddo
close(nut)
end subroutine
!**********************************************************
! ( 26 ) 
subroutine write_mloc_sub(mlocfile,incl,decl,geomag,aincl,adecl,idir,npt,oee,onn,oele,ai,ad)
character*(*) mlocfile
real incl,decl,geomag,aincl,adecl
integer idir
integer i,npt
real oee(npt),onn(npt),oele(npt),ai(npt),ad(npt)
integer nut
call open_or_not(nut)
open(nut,file=mlocfile)
  write(nut,*) npt
do i=1,npt
  write(nut,'(5f)')  oee(i),onn(i),oele(i),ai(i),ad(i)
enddo
close(nut)
end subroutine
!**********************************************************
! ( 27 ) 
subroutine write_mobs_sub(mobsfile,incl,decl,geomag,aincl,adecl,idir,npt,oee,onn,oele,ai,ad,mag,merr)
character*(*) mobsfile
real incl,decl,geomag,aincl,adecl
integer idir
integer i,npt
real oee(npt),onn(npt),oele(npt),ai(npt),ad(npt),mag(npt),merr(npt)
integer nut
call open_or_not(nut)
open(nut,file=mobsfile)
  write(nut,*) npt
do i=1,npt
  write(nut,'(7f)')  oee(i),onn(i),oele(i),ai(i),ad(i),mag(i),merr(i)
enddo
close(nut)
end subroutine

!**********************************************************
!********** this is the end of UBC-IO.f90 *****************
!**********************************************************